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Abstract 

We present an ab initio approach to solve the time-dependent Schrodinger equation to treat 
electron and photon impact multiple ionization of atoms or molecules. It combines the already 
known time scaled coordinate method with a new high order time propagator based on a predictor- 
corrector scheme. In order to exploit in an optimal way the main advantage of the time scaled 
coordinate method namely that the scaled wave packet stays confined and evolves smoothly towards 
a stationary state the modulus square of which being directly proportional to the electron energy 
spectra in each ionization channel, we show that the scaled bound states should be subtracted from 
the total scaled wave packet. In addition, our detailed investigations suggest that multi-resolution 
techniques like for instance, wavelets are the most appropriate ones to represent spatially the scaled 
wave packet. The approach is illustrated in the case of the interaction of an one-dimensional model 
atom as well as atomic hydrogen with a strong oscillating field. 
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I. INTRODUCTION 



During the last few years, substantial progress regarding the development of new XUV 
sources, has been made in two directions. On the one hand, high order harmonic generation 
has been used to produce attosecond pulses of which the duration is of the order of the 
characteristic time scale of the inner-shell electron dynamics in atoms and molecules [1]. On 
the other hand, free electron lasers [2] are now operating at unprecedentedly high peak inten- 
sities in the far-X-ray regime. These developments have opened the route to the exploration 
of non linear processes in the short-wavelength limit. At present, processes such as multi- 
photon multiple ionization of atoms and molecules are the focus of many experimental and 
theoretical studies with a view to understanding the subtle role of the electronic correlations. 

Within this context, there is clearly a need for reliable theoretical and numerical methods 
that provide accurate solutions of the time-dependent Schrodinger equation (TDSE). To 
this end however, it is necessary to overcome the following four main difficulties, (i) The 
continuum components of the wave packet expand in a rapidly increasing volume of space, 
thereby requiring very extended spatial grids or basis functions in order to avoid artificial 
reflections from the numerical boundaries, (ii) Increasingly large spatial phase gradients 
develop within the wave packet with time, demanding very dense grids or large basis sizes, 
(iii) Solving the TDSE on a spatial grid or in a basis of square integrable functions leads to 
a stiff system of equations which, in principle, makes explicit time propagators unstable. 
Finally, (iv) the direct extraction of the information on the multi-electron continua from 
the wave packet necessitates the knowledge of the asymptotic behavior of the corresponding 
wave function. 

The existing time dependent approaches have been mainly used to study single and 
double ionization of two-electron atoms and molecules by intense ultrashort radiation 
fields. In the low frequency regime where the calculations are extremely challenging, Smyth 
et al. [3] have developed a fully numerical method to solve the TDSE. It has provided 
valuable qualitative information on the role of the electronic correlations and the so-called 
rescattering process [4]. In the high frequency regime where one or two photons are 
involved in the ionization process, there are presently two types of treatment to solve the 
TDSE: the treatments based on standard methods of collision theory and the close-coupling 
approaches. In the former case, the wave packet is time propagated on an extended spatial 
grid during a period of time that is much larger than the pulse duration. The Fourier 
transform of the wave packet provides a scattered wave function which is then analyzed 
by means of time independent methods. Palacios et al. [5] use the Exterior Complex 
Scaling (ECS) technique which maps an outgoing wave into a vanishing wave outside 
a physically unaltered region allowing the extraction of the relevant information on the 
various ionization processes without the necessity of knowing the asymptotic behavior of 
the wave function associated to the multiple continua. Recently, Malegat et al. [6] applied 
the Hyperspherical R-Matrix with Semiclassical Outgoing Waves (HRM-SOW) method to 
calculate the various ionization yields. In this method, the scattered wave is propagated 
semiclassically with respect to the hyperradius, all the way to the asymptotic region where 
the various ionization channels are decoupled. 
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Many approaches based on a close-coupling method have been developed. They 
essentially differ by the way the information on the ionization processes is extracted from 
the wave packet. The most common way is to propagate the wave packet freely after the 
interaction with the pulse, until it reaches a spatial region where the ionization channels are 
assumed to be decoupled. It is then projected onto an uncorrelated product of Coulomb 
functions in each of the ionization channels [7-10]. Instead of using uncorrelated products 
of Coulomb functions, Ivanov and Kheifets [11] project the wave packet on continuum state 
wave functions obtained by means of the Convergent Close Coupling (CCC) method which 
takes into account electron correlations in an approximate way. A different procedure 
has been developed by Foumouo et al. [12]. Since the asymptotic behavior of the single 
continuum wave function is known, it is convenient to calculate the total probability 
for double ionization by subtracting the total probability for single ionization from the 
all-inclusive probability for breakup which in turn can be calculated without any reference 
to the boundary conditions. In order to calculate the total and partial probabilities for 
single ionization, they use the Jacobi-matrix method to generate a multichannel scattering 
wave function that describes the single continuum. The projection of the wave packet 
on this scattering wave function is performed just at the end of the interaction of the 
two-electron system with the radiation pulse. 

Finally, Lysaght et al. [13, 14] have recently initiated the development of a Time- 
Dependent i?-Matrix (TDRM) approach to describe complex multielectron atoms and 
atomic ions in intense ultrashort radiation pulses. This approach consists in time propagat- 
ing the atomic wave function in the presence of the radiation field both in the internal and 
external i?-matrix regions. 

The present approach combines the Time Scaled Coordinate (TSC) method with a 
high order fully implicit predictor-corrector scheme for the time propagation. The time 
dependent scaling of the radial electronic coordinates together with a phase transformation 
of the wave packet allow for "freezing" the spatial expansion of the wave packet in the 
new representation while removing fast oscillations due to the increasingly large spatial 
phase gradients that develop with time. This method is in fact equivalent to using a 
time-dependent basis that expands in the same way as the wave packet itself. This 
idea of time scaling the coordinates is not new and has been widely exploited in many 
different fields of physics. In 1979, Burgan et al. [15] studied the Schrodinger equation 
for a multidimensional quantum harmonic oscillator with time-dependent frequencies. By 
introducing an appropriate time- dependent scaling of the spatial coordinates, they were 
able to transform the problem to a free particle motion and to derive an exact analytical 
solution. Later on, Manfredi et al. [16, 17] introduced a time-dependent scaling of both 
space and time variables to "freeze" the expansion into a vacuum of both a one-dimensional, 
collisionless, two-species classical plasma and a quantum electron gas in planar geometry. 
In atomic and molecular physics, Solov'ev et al. [18] and later on, Ovchinnikov et al. [19] 
treated the Coulomb three-body problem, in particular ion-atom and atom-atom collisions, 
within a proper adiabatic representation by time scaling the internuclear distance. More 
recently, the TSC method has been used by Sidky et al. [20] and Derbov et al. [21] to treat 
the interaction of a model atom and molecule with an electromagnetic pulse and by Serov 
et al. [22, 23] to study electron impact single and double ionization of helium and more 
recently, double photoionization of two-electron atomic systems [24]. The TSC method is 
somehow an extension of a self-similarity analysis which has been introduced recently in 
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astrophysics [25]. By an appropriate scaling of all variables entering the equations governing 
the dynamics of a very large hydro dynamic system, the rescaled equations are identical to 
the original ones. This allows to define dual equivalent systems, the first one characterized 
by very long time and parsec length scales and the second one, characterized by very short 
time and small length scales allowing its study at the laboratory scale. Finally, the TSC 
method has been used to study the expansion of a Bose-Einstein condensate following the 
switch off of the trap [26, 27]. 

In the case of the interaction of an atom or a molecule with an electromagnetic pulse, 
the TSC method effectively confines the expansion of the scaled wave packet within a finite 
space of controllable size so that the evolution of this scaled wave packet can be followed over 
very long periods of time . Furthermore, it has been shown [21, 22, 28] that a long time after 
the end of the interaction of the atom or the molecule with the pulse, the energy spectrum 
of the ejected electrons is simply proportional to the modulus square of this scaled wave 
packet. The confinement of the scaled wave packet is due to three factors: the presence of an 
harmonic potential, the narrowing of the atomic potential and the increase of the effective 
mass of the electrons with time. This means that the effective de Broglie wavelength of 
these electrons decreases. In other words, the TSC method introduces different length 
scales in the problem. This has two important consequences. First, an optimal spatial 
description of the scaled wave packet requires multi-resolution techniques and second, it 
increases significantly the stiffness of the system of first order differential equations to solve 
for the time propagation of the wave packet. By this, it is meant that the time step rapidly 
decreases with increasing size of the system [29]. In this contribution, we describe a seventh 
order fully implicit predictor-corrector scheme. The predictor is the fifth order explicit 
method of Fatunla [30, 31] while the corrector is a seventh order fully implicit Radau 
method [34]. In principle, an implicit scheme requires solving large systems of algebraic 
equations at each time step. However, the accuracy of Fatunla's method is high enough 
[32, 33] to allow the use of an iterative procedure, the biconjugate gradient algorithm, 
to solve the large systems of algebraic equations at the corrector level. In other words, 
only matrix-vector products are needed, allowing a deep parallelization of the computer code. 

This contribution is organized as follows. In the first section after this introduction, we 
describe the TSC method in detail. For the sake of illustration, we consider the interaction 
of an one-dimensional system modeled by a Gaussian potential and interacting with a cosine 
square electromagnetic pulse. First, we examine the different reasons for the confinement of 
the scaled wave packet. Then, we study various spatial representations of this scaled wave 
packet and study its behaviour at various times after the pulse has ceased to interact with 
the model atom. In the next section, we describe our time propagation method. We first 
start with the explicit Fatunla's method and then discuss in detail the predictor-corrector 
scheme. The third section is devoted to the calculation of the energy spectrum. We derive 
an analytical expression in the case of our model atom and for atomic hydrogen. Results for 
both cases are presented and discussed in detail. The last section is devoted to conclusions 
and perspectives. Unless stated otherwise, atomic units are used throughout this paper. 



4 



II. THE TIME SCALED COORDINATE METHOD 



A. Outline of the method 



Our one-dimensional model that serves as an illustration for describing the TSC method 
consists of an electron initially bound in a Gaussian potential and interacting with a cosine 
square electromagnetic pulse. The TDSE that governs the dynamics of the electron is: 

i|U(x,t) = (H (x) + H I (x,t))^(x,t). (1) 
The atomic Hamiltonian H (x) is given by: 

H (x,t) = -~ + V(x), (2) 

where: 

V(x) = -V e-^ 2 . (3) 

By adjusting the parameters V and (3 that fix the depth and the width of the Gaussian 
potential, we can easily vary the number of bound states. In all the calculations we perform, 
we always assume that the model atom is initially in its ground state. Within the dipole 
approximation and in the velocity form, the interaction Hamiltonian Hj(x,t) writes: 

d 

Hj(x, t) = -iA f(t) sm(ut + <p) — . (4) 

A is the amplitude of the vector potential that is polarized along the x-axis. u is the 
frequency and <p the carrier phase. fit) is the pulse envelope defined as follows: 

1*1 <i 

/(*) = < (5) 

1*1 > 5 

The total pulse duration r = 2nn c /u where n c is an integer giving the number of optical 
cycles. The fact that n c is an integer is important since it ensures that the electric field has 
no static components. 

According to the TSC method [20], we introduce the scaled coordinate £ given by: 

(6) 




where R(cr) is an arbitrary scaling function. It is important to stress that a is just a 
parameter. In the following, we assume that it coincides with the time t. We write the 
scaled wave packet as follows: 

mt) = VRe^ 2 ^(R^t), (7) 



where the dot indicates the time derivative. The factor yR ensures that this scaled wave 
packet is correctly normalized. The phase transformation absorbs the fast oscillations of the 
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unsealed wave packet during its time evolution. The scaled wave packet $(£, t) satisfies the 
following TDSE: 



*«,«)• (8) 



Since the idea behind the TSC method is to build a time-dependent basis that expands 
in the same way as the wave packet, it is expected that if this expansion is accelerated, 
non inertial forces should appear. This explains the presence of the harmonic potential in 
the above TDSE. When this expansion occurs at constant velocity, i.e. when the scaling 
function is linear with time, the harmonic potential disappears. For R > 0, this potential 
confines the wave packet in a finite space. In fact, in the absence of external fields, the 
spectrum of the operator in square brackets in Eq. (8), becomes purely discrete [23]. An 
analysis of Eq. (8) shows that the scaling transformation (6) introduces an electron effective 
mass that is proportional to R 2 . Finally, we see that when R increases, the Gaussian 
potential narrows with time. This shrinking mainly affects the bound state components of 
the wave packet. Note that in the case of a Coulomb potential it is the effective electric 
charge that goes to zero with increasing values of R. 

Before analysing in more details the different factors that lead to the confinement of the 
scaled wave packet, let us examine the scaling function R(t). As stressed in [20], this function 
is arbitrary and chosen to facilitate the numerics. It must however satisfy a few constraints. 
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Figure 1: (Color online) Scaling function and pulse envelope as a function of time. For the scaling 
function, n = 3, t sc = initial and Roo = 0.025. The pulse has a sine square envelope and a total 
duration of 94 a.u.. 



It should be real, larger than one and equal to one from t = t initial corresponding to the 
beginning of the interaction until t = t sc where the scaling starts. In addition, for large time 
t, R should tend to R^t where R^ is what we call the asymptotic velocity. This ensures 
that for large times, the scaled wave packet becomes stationary. In the present calculations, 
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we define R(t) as follows: 

f 1, t < t sc 

R(t) = (9) 

[{l + [Roo(t-t sc )] n }n, t>t sc . 

This form with n = 4 has been used by several authors [20, 21, 28]. It leads to a smooth 
transition to the linear regime where the harmonic potential is switched off. The scaling 
function for n = 3, t sc = initial and Roo = 0.025 together with a sine square pulse envelope 
are shown as a function of time in Fig. 1. In the shaded region, between 10 and 60 a.u. of 
time, R(t) is significantly larger than zero. In that case, the confinement of the scaled wave 
packet results predominantly from the presence of the harmonic potential. For times t > 60, 
the harmonic potential is smoothly switched off but the confinement of the scaled wave 
packet subsists because of the electron effective mass which rapidly increases. As soon as 
the scaling starts (at t = t sc ), the atomic bound state wave functions start to shrink because 
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Figure 2: Effective potential namely the sum of the scaled Gaussian potential and the harmonic 
one together with the scaled ground state wave function versus the scaled variable £. Three times 
are considered: (a) t = t sc corresponding to the time where scaling is switched on (b) t = t sc + 15 
a.u. and t = t sc + 80 a.u.. The asymptotic velocity is equal to 0.1 a.u. and n = 4 (see Eq. (9)). 

of the narrowing of the atomic potential. In the case of the Gaussian potential it is easy to 
show that the width is inversely proportional to R. This is illustrated in Fig. 2 where we 
show the scaled ground state wave function as well as the effective potential which is the sum 
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of the scaled Gaussian potential and the harmonic one as a function of the scaled variable 
£, for three different values of t: t sc , t sc + 15 a.u. and t sc + 80 a.u.. Note that at t — t sc , the 
effective potential reduces to the Gaussian potential. The fact that the scaled ground state 
wave function shrinks rapidly is due to the relatively high value of the asymptotic velocity 
which, in the present case, is equal to 0.1 a.u.. As n = 4 (see Eq. (9)), t — t sc + 80 a.u. 
corresponds to a time for which the scaling function has reached its linear regime. It is 
important to mention that the results given in Fig. 2 require only a partial diagonalization 
of the atomic Hamiltonian and no time propagation. Therefore, they provide an easy way 
to control the fineness of the spatial discretization necessary to maintain the accuracy of the 
time propagation. 



B. Spatial representation of the scaled wave packet 

The optimal way of describing the wave packet in space is based on a multi-resolution 
analysis [35]. The general idea is to define different resolution levels in various regions of 
space through the introduction of several grids with a density of mesh points that increases 
from one grid to the next one in the spatial regions of interest. These techniques will 
be analysed in detail in a forthcoming publication. Here, we use two different spectral 
methods. The first one consists in developing the wave packet on C 2 integrable functions 
namely Hermite-Sturmian functions in the case of our one-dimensional model and Coulomb- 
Sturmian functions in the case of atomic Hydrogen. The second method uses a basis of 
B-splines built on a non-uniform grid with an exponential sequence of breakpoints. These 
two methods that are far from being optimal, have the merit to be easily implemented and 
to clearly show that one cannot dissociate the spatial discretization problem from the time 
propagation. 

Let us now briefly describe our first spectral method. In the case of our one-dimensional 
model, we expand the total wave packet in a finite basis of Hermite-Sturmian functions as 
follows: 

N 

9(x,t) = Y,an(t)(ft(x), (10) 

n=0 

where a n (t) is the expansion coefficient and <p%(x), the Hermite-Sturmian function given by: 

_ i 

<fi(x)=(2 n n\l?) \~^H n {^x). (11) 



The elements of the matrices associated to all operators present in the scaled and unsealed 
Hamiltonians can be calculated analytically except for the Gaussian potential. In this 
latter Gauss-Hermite quadrature provides exact results for a sufficient number 

of abscissae, a is a dilation parameter that determines the resolution of the basis. 
Indeed, a large value of a gives a good description of a wave packet which exhibits sharp 
variations close to the origin. In that case however, a large value of N, the number 
of basis functions, is necessary if the extent of this wave packet is significant. By con- 
trast, a small value of a allows a good description of the wave packet over much larger 
distances but N has to be extremely large again if sharp variations of the wave packet occur. 
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In the case of atomic hydrogen where we use spherical coordinates, we write the wave 
packet as follows: 

*(?,t) = mW^vm), (12) 

n,l,m 

where Yi !Tn (6,4>) is a spherical harmonic and S^^r) the Coulomb- Sturmian function given 
by: 

SUr) = r' +1 e— L?+ 1 _ 1 (2«r). (13) 

N£ t is a normalization factor (see [12] for details) and L^^_ 1 (2/cr) a Laguerre polynomial. 
The index n varies from / + 1 until oo. The Coulomb Sturmian functions form a complete 
and discrete set of C 2 integrable functions that are the exact solutions of the stationary 
Schrodinger equation for a single electron in the Coulomb field of a nucleus of charge 
Z for selected values of k. As a result, these functions are well adapted to describe the 
energy spectrum of atomic hydrogen as well as its behavior in presence of an external 
field. In fact, the matrices associated to the corresponding Hamiltonian are banded with 
a narrow bandwidth (three diagonals). As in the case of the Hermite Sturmian functions, 
a given basis of Coulomb Sturmian functions is characterized by a fixed value of the 
dilation parameter k. An interesting idea in the spirit of the multi-resolution approaches 
is to consider a set of different values of k within a given basis in order to take into 
account the various length scales in the problem. Despite the fact that the introduction 
of different values of k makes the basis numerically overcomplete thereby requiring the 
elimination of the linearly dependent eigenvectors of the overlap matrix, this idea turned 
out to be extremely successful to generate the singly and doubly excited states of helium 
[12, 36-39]. Note that in the case of atomic hydrogen, the time dependent scaling of the 
radial coordinate is equivalent to introducing a time dependent k. 

In our second method to treat our one-dimensional model, we expand the wave packet in 
a basis of B-splines [40] 

N 

i=i 

where Bf(x) is a B-spline of order k. In the present calculations, k = 7. In order to calculate 
the B-splines, we use an exponential sequence of breakpoints. In practice, we proceed as 
follows. We adjust the asymptotic velocity and the time t sc at which the scaling starts so 
that the scaled wave packet is confined in a relatively small interval [-x max , +x max ] . We 
then define two symmetrical exponential sequences of breakpoints in the intervals [-x ma _ x , 0] 
and [0 , +x max ]. In the interval [0 , +x max ] for instance, the breakpoints £j are given by: 



i-l 



•'•u,ax I e7-1 j •' I V. (15) 

Typically, we have x max = 35, 7 = 5 and N of the order of 100. In this case, breakpoints 
accumulate symmetrically around in order to describe properly the shrinking of the 
bound states. However, this accumulation of breakpoints introduces high frequencies in the 
problem. In this B-spline basis, the energy spectrum of the unsealed atomic Hamiltonian 
contains very high positive eigenenergies of the order of 3000 a.u.. If the wave packet is 
time propagated in the B-spline basis, each component of the B-splines contains these high 
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frequencies making the problem extremely stiff and leading to a dramatic decrease of the 
time step. Note that this difficulty can be avoided by propagating the scaled wave packet 
in the atomic basis i. e. the basis in which the unsealed atomic Hamiltonian is diagonal. In 
that case however, it is necessary to solve a generalized eigenvalue problem that could be 
time consuming in the case of a more complex system like He or H2. 



C. Evolution of the scaled wave packet 

Before examining the time evolution of scaled wave packets, it is instructive to analyze 
the behavior of the unsealed ones. We first consider the case of a Gaussian potential with 
Vq = 1 a.u. and (5 = 1 a.u.. This potential has only one bound state the energy of which 
is equal to -0.477 a.u.. This model atom interacts with a cosine square laser pulse of peak 
intensity J pcak = 10 13 W/cm 2 and frequency u = 0.7 a.u.. The total duration of the pulse is 
6 optical cycles. In Fig. 3, we show the real part of the wave packet at time t = tfi na i i.e. 
at the end of the pulse (dark grey curve) and 300 a.u. of time later (light grey curve). We 
actually represent the ionized wave packet i. e. the total wave packet without its bound state 
component. We clearly see that at time t = tfi na i + 300 a.u., the real part of the ionized wave 
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Figure 3: Real part of the unsealed ionized wave packet resulting from the interaction of our one- 
dimensional model atom with a cosine square pulse of 10 13 Watt/cm 2 peak intensity and 0.7 a.u. 
photon energy. The total duration of the pulse is 6 optical cycles. The Gaussian potential depth 
Vq = 1 a.u. and (3 = 1 a.u.. The real part of the ionized wave packet is shown at time t = tfi na i i.e. 
at the end of the pulse (dark grey curve) and at time t = tfi na i + 300 a.u. (light grey curve) . 

packet exhibits a much larger number of oscillations than at t = t &na i, the end of the pulse. 



10 



In particular, we clearly see the presence of a chirp for t = £fi na i + 300. As stressed by Sidky 
et al. [20], this results from the phase gradients that rapidly develop over large distances or 





Figure 4: Real part of both the scaled (dark grey curve) and unsealed (light grey curve) ionized 
wave packets resulting from the interaction of our one-dimensional model atom with the same pulse 
as in Fig. 3. The real part of these wave packets is calculated at time t = ifi na i + 250 a.u. and is 
represented as a function of the position: x in the case of the unsealed wave packet and £ in the 
case of the scaled wave packet. The parameters of the Gaussian potential are the same as in Fig. 
3. In the case of the scaled wave packet, the asymptotic velocity Roq = 0.01 a.u., the parameter n 
of the scaling function (see Eq. 9) is equal to 4 and the scaling is switched on (a) at the end of the 
interaction and (b) at the beginning of the interaction of our model atom with the pulse. 



in other words, from the fact that the front edge of the wave packet is moving faster than 
the inner part. Let ^f(x, £fi na i) be the wave packet created at the end of the pulse. At any 
later time t, we have: 

/oo 
dx' 0(x,x',t)* (affinal), (16) 
-oo 
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where the Green function Q(x,x',t) for a free electron is given by: 



Q(x,x',t) = ^ t e^ 2 ^. (17) 

This means that the phase increases quadratically with the distance. In the present calcula- 
tions, we need to use a basis of 1000 Hermite Sturmian functions of parameter a = 0.01 to 
accurately reproduce all the oscillations over about 300 a.u. around the origin. In Fig.4a, 
we compare for t = tg na i + 250 a.u., the real part of both the scaled and unsealed ionized 
wave packets for the same pulse and Gaussian potential parameters as in Fig. 3. In this 
particular case, the scaled ionized wave packet is obtained as follows. The unsealed wave 
packet is first propagated until t = tfi na i, the end of the interaction of our model atom with 
the pulse. At time t = tfi na i, the bound state contribution is subtracted from the total wave 
packet and scaling is switched on. At time t = tfinai + 250 a.u., we clearly see on Fig. 4a that 
the scaled wave packet represented as a function of the position £ is confined compared to 
the unsealed one. In addition, the number of oscillations is already reduced. As discussed 
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Figure 5: Real part of the scaled wave packet resulting from the interaction of our one-dimensional 
model atom with a cosine square pulse of 5 x 10 14 Watt/cm 2 peak intensity and 0.3 a.u. photon 
energy. The total pulse duration is 6 optical cycles. The parameters of the Gaussian potential 
are the same as in Fig. 3. The time dependent scaling is switched on at the beginning of the 
interaction with the pulse. The asymptotic velocity Roo is equal to 0.08 a.u. while the parameter 
n of the scaling function (see Eq. 9) is equal to 4. The real part of the scaled wave packet is 
represented as a function of £ at (a) t = tg na i and (b) t = tfi na i + 1930 a.u.. 

above, the quadratic increase of the phase with the distance is canceled by the phase trans- 
formation (7) of the wave packet. In principle, the time scaling of the coordinates may start 
at any time. In Fig. 4b, we consider the same case as in Fig. 4a but the time scaling is now 
switched on right at the beginning of the interaction with the pulse. We clearly see that the 
confinement is slightly stronger and that the number of oscillations is significantly reduced. 
Indeed, beyond ±120 a.u., all the oscillations of weak amplitude present in Fig. 4a have 
disappeared. The comparison of Figs 4a and 4b also shows that the scaled wave packet has 
not yet reached a stationary state. In Fig. 5, we show the scaled wave packet resulting from 
the interaction of the same one- dimensional model atom as before with a pulse of 5 x 10 14 
Watt/cm 2 peak intensity and 0.3 a.u. photon energy. The total duration of the pulse is 6 
optical cycles. In this particular case, we use B-splines to describe the scaled wave packet 
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with two different exponential sequences of break points. The first grid with 7 = 5 (see 
Eq. 15) and 200 B-splines is used during the interaction of the model atom with the pulse. 
The second grid with 7 = 5 and 400 B-splines is used after the interaction with the pulse 
i.e. for t > ifinai- In Fig. 5, we show the real part of the scaled wave packet at the end 
of the pulse (Fig. 5a) and 1930 a.u. of time later when the scaled wave packet becomes 
quasi stationary. The fast oscillations have now been completely removed. It is important 
to stress that although the number of B-splines used is rather small, both grids allow us 
to describe accurately the shrinking of the ground state. In fact, for time t ^> ifi na i when 
the scaling function becomes linear with time, the scaled TDSE is identical to the original 
TDSE within a simple linear scaling of the spatial coordinates. This clearly shows that a 
multi-resolution approach in which the density of break points increases linearly with time 
around the origin makes sense. Note that in the present case, the accumulation of break 
points around the origin (right from the beginning of the propagation) increases significantly 
the stiffness of the system, a problem which should be avoided by means of multi-resolution 
techniques. The stiffness problem is discussed in the next section. 

III. TIME PROPAGATION 

Solving accurately the TDSE usually requires the representation of the solution on 
large or/and dense grids or in large bases. In all cases, we deal with large systems of 
coupled first order differential equations which are well known to be stiff [41]. This means 
that the step size decreases as the dimension of the system increases. The origin of the 
stiffness is clear: by increasing the size or the density of the grid or the size of the basis, 
the Hamiltonian generates large positive eigenenergies, which are responsible for strong 
oscillations in the solution of the TDSE. It is thus the largest positive eigenvalue which 
controls the step size. In fact, the stiffness of the system may lead to the instability of the 
time propagation scheme as well as to inaccurate high energy components of the solution 
[42]. Two approaches frequently used to overcome this problem are implicit schemes for 
solving the TDSE and the propagation of the TDSE in the atomic basis and possibly, within 
the interaction picture. The first method requires typically the solution of large systems 
of linear equations at each integration step. It is however important to stress that implicit 
schemes actually solve the stability problem but not necessarily the inaccuracy problem for 
the high energy components of the solution. In the second method, the time integration is 
achieved by means of explicit algorithms, which have been proved to be very stable for the 
solution of the TDSE in the atomic basis where the atomic Hamiltonian is diagonal. These 
algorithms only need matrix-vector products. However, the representation in the atomic 
basis requires the full diagonalization of the Hamiltonian before starting the integration. In 
any case, the computational cost increases dramatically with the size of the system. 

It is therefore desirable to have an explicit algorithm suitable for the direct solution of stiff 
TDSE. Such a method does exist and was proposed more than thirty years ago by Fatunla 
[30, 31]. It has been successfully implemented for the description of single ionization of 
atoms by strong oscillating fields [32, 33]. In this method that takes into account the 
intrinsic frequencies of the system, the wave function is expressed in terms of oscillating 
functions. This leads to a simple recursive formula for the time propagation with a controlled 
error. At each integration step, only matrix-vector products are therefore needed. In the 
two following subsections, we describe the most important features of Fatunla's algorithm 
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and its implementation and show how its accuracy can be significantly improved within a 
predictor-corrector scheme. 



A. Fatunla's explicit scheme 

We start with the general matrix form of the TDSE using a spectral or a grid represen- 
tation, 

iB^=H(f)*, (18) 

with H(£) the matrix representation of the Hamiltonian, and B the overlap matrix, which 
is the identity in a grid representation or in an orthonormal basis. Truncation of the basis 
or of the grid leads to a m-dimensional first order differential equation, 

* = /(<,*), (19) 
where f(t, \I>) = — iB 1 H(t)* is in general a complex m-dimensional function. 

The stiffness of equation (19) leads to a solution &(t) which is an oscillating function. In 
a given interval [t n ,t n+ \], t n+ i — t n + h, with h a small number, &(t) is approximated by 
the function 

F(t) = (I - e Qlt )a - (I - e-^^b + c, (20) 

with I the identity matrix, f2j = diag(a>f\ . . . , oom), i = 1,2, and a, b, c constant vectors. 
The complex numbers oo[ l \ . . . , oom , i — 1,2 are called stiffness parameters. Assuming that 
F(t) coincides with &(t) at t n and t n+ i, that F'(t) coincides with f(t,&) at t n , and that 
F"(t) coincides with f'(t,*&) at t n , the solution Vl/ n +i = &(t n+ i) at t n+ i can be expressed 
recursively in terms of = *(£ n ), f n = f(t n , * n ) and f^> = df/dt\ t=tn according to 

* n+1 = * n + R/ n + S/«. (21) 

R and S are diagonal matrices which can be written in terms of the stiffness parameters: 

R = ft 2 $-fi 1 3, S = $ + H, (22) 

where <fr and S are diagonal matrices of which the nonzero entries are 

e J1)h - 1 

^ = (1)771)7 W, ^ 
uj) '{uj) + toy) 



and 

(24) 



Notice that if a stiffness parameter, uf^ {oof"*), vanishes the associated matrix element reads 

*.=4> ( 3 -=-4>V < 25 > 
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The recursive relation (21) depends on the so far unknown stiffness parameters. These can 
be written in terms of the function f(t n , & n ) and its time derivatives k = 0, 1, 2, 3, at 

tn [31], 



4 2) = u,® + A, (26) 

where and E i} i = 1, . . . , m are given in terms of the respective components of 
fc = 0,l,2,3atf = t n by 

AO) ,(3) _ ,(i) A2) 

/-) J in J in J in J in ■ i (07\ 



J in J in J in J \ 



in 



M) A3) _ ,(2) ,(2) 

rp J in J in J in J in ■ i (OQ'\ 

^ ~ ,(1) ,(1) _ ,(0) ,(2) > I- 1, 
J in J in J in J in 

provided that the denominator of the previous expressions is nonzero. It is important to 
note that in the present case, the time-dependent scaling function (9) is considered as a 
parameter. As a result, the successive time derivatives of this function are not taken into 
account in the calculation of _Dj and E^. 

The ith component of the local truncation error at t — t n+ i, that is the difference between 
the exact solution at i n+1 and the numerical solution, is given by (we ignore here the index 
i of the components for the sake of clarity) [31]: 

o! uj\ + U2 

-( u * lUj2 + UJlUJ *)f(°)] + 0(h 6 ) 

-u^ul - u x u 2 + u 2 2 )f^] + 0(h 6 ). (29) 

The implementation of the recursion (21) is now rather simple. It requires the calculation 
of the function f n and its derivative at each value of t n . For the stiffness matrices 
fli and f2 2 , and thus also for the matrices R and S, the derivatives and fffl are also 
needed. f2i and f2 2 have to be calculated in principle at each integration step, since they 
characterize the local frequencies of the solution *&(t). The truncation error (29) can be 
used to control the size of the integration step, e.g., by imposing a boundary criterion for 
|T n |. For this also, the derivative fffl must be provided. 



The stiffness parameters carry the intrinsic information of the natural oscillations of the 
system. Therefore, the time step is expected to be rather large compared with standard 
explicit methods, like Runge-Kutta. In order to illustrate this, we solve without any time- 
dependent scaling of the coordinates, the TDSE (18) for our one-dimensional model atom 
(Vo = 1 and f3 — 1) interacting with a 6-cycles cosine square pulse of 10 13 Watt/cm 2 peak 
intensity and 0.7 a.u. photon energy (same case as in Fig. 3). The time propagation is 
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Figure 6: Comparison of the time steps used in Fatunla's method (dashed line) and a method 
based on a fifth order embedded Runge-Kutta formula (dotted line) to solve the TDSE for the 
same case as in Fig. 3. The time propagation is performed in a basis of 500 Hermite Sturmian 
functions of dilation parameter a = 0.5 a.u.. For the sake of completeness, we also show the pulse 
envelope. 

performed in the Hermite Sturmian basis by means of two approaches : Fatunla's method 
and a fifth-order embedded Runge-Kutta formula [43] . The Hermite Sturmian basis contains 
500 functions with the dilation parameter a = 0.5 a.u.. In Fig. 6, we compare the time 
variations of the time step which is automatically adjusted in both methods. In the case 
of Fatunla's method, this is done according to the condition 10~ 16 < T n < 10~ 12 for the 
truncation error. In the case of the embedded Runge-Kutta formula, the time step is adjusted 
by comparing the results obtained with the fifth order formula and with a fourth order one 
which uses the same intermediate mesh points. We clearly see that the time step is bigger 
in the case of Fatunla's method when compared to the Runge-Kutta embedded formula. 
However, the relative error on the conservation of the norm of the wave packet is, in this 
particular case, two orders of magnitude better in the case of the Runge-Kutta embedded 
formula. In fact, by increasing the value of the dilation parameter a as well as the number 
of basis functions, the stiffness of the system of equations becomes more pronounced. In 
these conditions, Fatunla's method becomes much more efficient in terms of the magnitude 
of the time step and the relative error on the conservation of the norm is of the same order 
of the one obtained with the Runge-Kutta embedded formula. 

Fatunla's method has also been tested in much more demanding cases namely the 
interaction of atomic hydrogen with intense low frequency pulses. The corresponding TDSE 
has been solved in a Coulomb Sturmian basis. The results obtained with Fatunla's method 
are in very good agreement with those obtained with a fully implicit Radau method of 
order 7 (see reference [33, 41] for details). In fact, Fatunla's method turned out to be 10 
times faster than the fully implicit method. However, the relative error on the conservation 
of the norm of the solution is hardly lower than 10~ 5 by contrast with the fully implicit 
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Figure 7: Comparison of the time steps used in Fatunla's method to solve the scaled TDSE (light 
dotted line) and the unsealed TDSE (thick dotted line) in the case of the interaction of our one- 
dimensional model atom (Vo = 1 and (3 = 1) with a cosine square pulse of 10 13 Watt/cm 2 peak 
intensity and 0.7 a.u. photon energy. The total pulse duration is 40 cycles. The scaled and 
unsealed TDSE are solved in a Hermite Sturmian basis. In the case of the scaled TDSE, 400 
Hermite Sturmian functions are used. These functions have a dilation parameter a = 0.05 a.u.. 
The parameter n of the scaling function is equal to 4 and the asymptotic velocity i?oo = 0.01 a.u.. 
For the unsealed TDSE, a basis of 2000 functions with a dilation parameter a = 0.008 a.u. is 
necessary. For the sake of completeness, we also show the pulse envelope. 

method where the norm is perfectly conserved. 

Let us now examine how Fatunla's method performs in the case of the scaled TDSE. We 
consider the interaction of our one-dimensional model atom (Vo = 1 a.u. and ft — 1 a.u.) 
with a cosine square pulse of 10 13 Watt/cm 2 peak intensity and 0.7 a.u. photon energy. 
The total pulse duration is equal to 40 optical cycles. In Fig. 7, we show how the time step 
used by Fatunla's method varies with time for the scaled and unsealed TDSE. We work in 
the Hermite Sturmian basis. In the case of the scaled TDSE, 400 basis functions with a 
dilation parameter a = 0.05 a.u. are sufficient. The parameter n of the scaling function 
is equal to 4 and the asymptotic velocity Roo = 0.01 a.u.. It is important to note that 
in this case, it is necessary to integrate over about 4000 a.u. of time to reach the region 
where the scaled wave packet becomes static. By contrast to the unsealed equation, we 
clearly see that the time step needed for the time propagation of the scaled wave packet 
increases significantly as soon as the interaction with the pulse has ceased. In fact, around 
t = 4000 a.u., the time step has a value around 10 a.u.. This results from the subtraction 
of the scaled bound state at the end of the pulse and from the slow time evolution of the 
ionized wave packet in the scaled representation. It is important to stress here that the 
shrinking of the bound states will never stop whereas the ionized part of the wave packet 
becomes static at large times. It is therefore important to subtract the scaled bound states 
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from the total wave packet. Note that the subtraction has to be performed when both the 
interaction with the pulse and the harmonic potential have disappeared. 

From the previous discussion, it turns out that Fatunla's method with adaptive step size 
is particularly well adapted to the solution of the scaled TDSE. In very stiff problems, the 
time step is much bigger than in the case of usual explicit Runge-Kutta methods. However, 
in all cases we have treated so far, the relative error on the conservation of the norm is of 
the order of 10~ 5 even when we force the time step to be very small. In many cases, this 
accuracy is sufficient but, there are always cases where a higher accuracy is needed. In the 
next subsection, we show that this accuracy problem is solved by using a predictor-corrector 
scheme in which Fatunla's method is the predictor. 



B. Predictor-corrector scheme 

Predictor-corrector (P-C) methods are pairs of an explicit and an implicit multistep 
method where the explicit formula is used to predict the next approximation and the implicit 
formula to correct it. The order of the implicit method is usually the same or higher than 
the order of the explicit method. In the present case, the predictor is Fatunla's method 
which is of order 5. The corrector is a fully implicit 4-stage Radau method of order 7 [34]. 
The implementation of this implicit method within the P-C scheme follows Refs. [44] and 
[45]. It is based on diagonally implicit iterations that have two main advantages: it preserves 
the favourable stability characteristics of the fully implicit Radau method and it is suitable 
for use on parallel processors. In the following, we give a brief outline of the method. Let 
t n and t n+ i, be two consecutive times at which we want to calculate the wave packet. The 
solution *& n+ i = *&(t n+ i) of the TDSE is obtained from \I> n through the following relation: 

4 

* n+ i = *„ + /i^6 J H(t n + c l /i)Y J , (30) 
i=i 

where h = t n+ i — t n and 6« and q are coefficients that define Radau's method, ti = t n + C{h 
with % varying from 1 to 4 are intermediate times in the interval [t„,t„+i] with t 4 = t n+ \ 
(04 = 1). In the case of an orthonormal basis, H is the Hamiltonian matrix. When the 
basis is not orthonormal, H represents the Hamiltonian matrix multiplied on the left by the 
inverse of the overlap matrix. Yj is an estimation of *f?(ti) which, within the fully implicit 
4-stage Radau method, is obtained by solving the following system: 













H 




H 











fluH(ti) ••• a 14 H(t 4 ) \ /Yx 



a 4 iH(ti) • • • a 44 H(t 4 ) 









!v,) 



(31) 



where the coefficients are given. They define like bi and q the implicit 4-stage Radau 
formula. If H is a m x m matrix, the dimension of the above system is 4m. Usually, when 
an implicit method is the corrector in a P-C scheme, the vector that contains the Yj in the 
right hand side of Eq. (31) is provided by the predictor and system (31) is solved iteratively. 
Those iterations are called explicit since they only require matrix-vector multiplications. 
However, for stiff problems, such an explicit iterative process does not always converge. 
Instead, we introduce the diagonal matrix D=diag(c?n, d 22 , rf 33 , <i 44 ) and rewrite system (31) 
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as follows: 



/ Y? \ 


/dnH(*i) 








V o 



















/ Y? 



■ ■ (i 44 H(t 4 ) 
(a n - dn)H(t!) 

a 4 iH(ti) 



ai 4 H(i 4 ) 
(a 44 - d 44 )H(i 4 ) 



\Yrv 



,(32) 



where the superscript j gives the order of the iteration, j — corresponds to the result 
obtained with Fatunla's method. Although systems (31) and (32) are perfectly equivalent 
since matrix D is subtracted on both sides of Eq. (31), we have now to solve four m x m 
systems of equations at each iteration. Those iterations are called implicit. Note that due 
to the diagonal structure of matrix D, the four systems can be solved in parallel. Such an 
implicit iterative process converges rapidly requiring a few iterations. The choice of matrix 
D is in principle arbitrary. In the present case, it is calculated in order to preserve the 
stability properties of the fully implicit Radau method. See Ref. [44] for more details. The 
present P-C method allows a perfect conservation of the norm. Furthermore, we checked 
that the four systems of equations can be easily solved by means of the biconjugate gradi- 
ent algorithm which only requires matrix vector multiplications. Therefore, although the 
present P-C algorithm is of an implicit nature, it is easy to implement on parallel processors. 



IV. ELECTRON ENERGY SPECTRUM 

One of the main advantages of the TSC method is the fact that the electron energy 
spectrum may be expressed directly in term of the scaled wave packet at large times. Here, 
we show it explicitly in the case of our one-dimensional model and for atomic hydrogen. 



A. Analytical expression for the electron energy spectra 

1. One- dimensional model 

The unsealed wave packet that is the solution of the TDSE (1) can always be written as 
follows: 

/OO 2 
c^M^e-^dk, (33) 
■°° 

where if) n (x) is the wave function associated to a bound state of our model atom with E n its 
energy. a n (t) represents the probability amplitude for the system to be in the corresponding 
bound state. ipk{x) is the wave function associated to a continuum state of energy E = k 2 /2 
while Ck(t) is the corresponding probability amplitude, k is the wave vector: a positive value 
of k corresponds to a wave propagating to the right while a negative value corresponds to a 
propagation to the left. The continuum wave functions satisfy the following orthogonality 
relation: 

/oo 
r k (x)^(x)dx = 5(k-k'). (34) 
-oo 
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In order to calculate the electron energy spectrum given by \ck(t — > oo)| 2 , we let x — > ±00 
in Eq. (33). In this limit the first term of the right hand side of this equation goes to zero 
and, 

~> 00) = A^f kx - (35) 



Note that the limit x — > —00 leads to the same expression since k is replaced by —k. As a 
result we have 

V(x ->■ ±00, t) = —= / Cki^e^'^Jdk. (36) 

V27T J-oo 

In order to calculate the spectrum, let us take the limit of the above integral for t — > 00. By 
using the stationary phase theorem, we obtain: 

1 

V(x -» ±00, t ->• 00) = c fc (i ->■ oo)^=e 1 2i. (37) 

Vit 

where = x/t is the stationary phase point. From Eqs. (6) and (9), we write 

k = j = R 0O (l- t -f)t. (38) 
Using Eq. (7) that relates the scaled wave packet to the unsealed one, we finally get: 

c k (t ^oo) = Ji( — + ^)$( — ^_ q , t^oo], (39) 

which establishes a direct link between the probability amplitude for the electron to be in 
the continuum with an energy k 2 /2 and the scaled wave packet at large times. Below, we 
discuss results for the electron energy spectrum and compare them with those obtained 
by solving the unsealed TDSE. In this latter case, we calculate the energy spectrum by 
projecting the final wave packet on continuum pseudostates obtained by diagonalizing the 
atomic Hamiltonian. Note that we get the same result by projecting onto plane waves 
because the phase shift introduced by the short range Gaussian potential is negligible. 



2. Atomic hydrogen 

As in the previous case, the unsealed wave packet can be written as follows: 

poo 2 

tf(f,f) = V M^n.iMJWMe-^' + V / c lk (t)R l (kr)Y l>m (0,<l))e' i ^ t kdk, (40) 

n,l,m l,m 

where ip n ,i{i r ) is the radial wave function of an hydrogen bound state of energy E n and a n j 
the corresponding probability amplitude. Y^ m (9, </>) is a spherical harmonic that is function 
of the angular coordinates of vector r. Ri(kr) is the radial regular Coulomb function where 
k is the magnitude of the electron momentum [46]. Ci k {t) is the probability amplitude for 
the electron to be in the continuum. The asymptotic form of the radial regular Coulomb 
function is 

Ri{kr) ->■ NAk) sin ( kr - -In + ^ln(2A;r) + aA . (41) 
r^oo fcr V 2 A; / 
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Ni(k) is a normalization factor. In the present case, it is equal to 2\fk because Ri(kr) is 
normalized per unit energy. 07 is the Coulomb phase shift given by 



<r, = argr(Z + l--), 



(42) 



where T(x) is the gamma function. In the limit of large distances, the unsealed wave packet 
reduces to the following expression: 



r — ^no * * 



Lm 



cikit) e 1 2 



3 i(fcr-^+<Ti + iln(2AT)) _ -i(kr- l -f+a l +\\n{2kr)) 



dk. 



Note that the sine function in the asymptotic expression (41) of the regular Coulomb func- 
tion has been replaced by a sum of two complex exponentials that describe pure ingoing 
and outgoing spherical waves at large distances. In order to calculate the electron energy 
spectrum, we now examine the limit of the previous expression for t — > 00. Before applying 
the stationary phase theorem, let us rewrite expression (43) as follows: 

^ c lk (t) ei^'f) e tt K + x+^<2*r)) dL (44) 



r — ^no * * 



Lm 



ir 



In fact, it is easy to show that the ingoing spherical wave present in expression (43) does not 
contribute to the integral in the limit t — > 00. In integral (44), the stationary phase point 
k = k is the solution of the following equation: 

1 - ln(2£x) 



-k + 



0. 



(45) 



t kH 

To a good approximation we can replace k by r/t in the third term of the left hand side of 
the above equation. As a result we obtain: 



k = 



ln(2V) 



t 



(46) 



If k e represents the velocity of the electron at large distances, we can show by means of 
classical mechanics, that for large times, 



(47) 



Note that for large times, k , the stationary phase point coincides with k e . By using the 
stationary phase theorem and Eqs. (6), (7) and (9), we obtain: 



c lko (t ->■ 00) 



+ 



^sc \ 

W)J 



G) 



e -i(<r,-ir/2) 



cxp 



'ln(2r 2 /t) + 2^) 



+ 



ln(2r 2 /f) 
r 2 /t 



X $; 



-, t — > OO 



R„(t-t K r -j- (48) 

which for a given value of the angular momentum /, establishes a link between the proba- 
bility amplitude for the electron to be in the continuum with an energy E = A^/2 w k%/2 
and the corresponding /-component of the scaled wave packet and hence the electron energy 
spectrum. Our results for the energy spectrum are compared with those obtained by 
projecting the unsealed wave packet on Coulomb functions. 
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B. Results and discussion 



As a proof-of-principle, we show in this section, that the present method provides very 
accurate results for the electron energy spectra at the expense of less computer resources 
than with the same propagation method without scaling. Here, we calculate two different 
electron energy spectra in rather demanding physical situations. First, we consider the 
ionization of our one-dimensional model atom with an intense low-frequency field. The 
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Figure 8: (Color online) Electron energy spectrum resulting from the interaction of our model atom 
with a strong low frequency pulse. The Gaussian potential parameters Vq = 4 a.u. and f3 = 0.1 a.u. 
are such that this potential supports seven bound states. The cosine square pulse of peak intensity, 
10 16 Watt/cm 2 and frequency u = 0.5 a.u. has a total duration of 8 optical cycles. The full blue 
line is the result obtained by solving the unsealed TDSE with a basis of 3000 Hermite Sturmian 
functions and a dilation parameter a = 0.01 a.u.. The red full line is the result obtained by solving 
the scaled TDSE with a basis of 1500 Hermite Sturmian functions and a dilation parameter a = 0.1 
a.u.. 

Gaussian potential parameters (Vo = 4 a.u., and (3 = 0.1 a.u.) are such that it supports 
seven bound states, the ground state energy being equal to -3.572 a.u.. The cosine square 
pulse has a total duration of 8 optical cycles with a peak intensity of 10 16 Watt/cm 2 and a 
frequency of 0.5 a.u.. The time propagation of both the unsealed and scaled wave packets 
has been performed by means of the predictor-corrector method described above. The 
unsealed TDSE has been solved by expanding the wave packet in a basis of 3000 Hermite 
Sturmian functions with a dilation parameter a = 0.01 a.u.. The result (the full blue line) 
is shown in Fig. 8. We clearly see that the energy spectrum becomes noisy already around 
an electron energy of 3 a.u.. The reason is the following. Since fast electrons are emitted, 
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the problem of the reflection of the unsealed wave packet by the numerical boundaries is 
a crucial issue. In order to overcome this problem, it is necessary to use a small value 
of the dilation parameter. In that case however, the density of positive energy states we 
obtain by diagonalizing the atomic Hamiltonian gets very small at high electron energy, 
thereby leading to the noisy behavior of the energy spectrum above 3 a.u.. The scaled 
TDSE has been solved by using a smaller basis of 1500 Hermite Sturmian functions with 
a dilation parameter of a = 0.1 a.u.. The correct spectrum (red full line) is obtained over 
more than 10 orders of magnitude after having time propagated the scaled wave packet 
over 5000 a.u. of time. The origin of the discrepancies observed at low electron energies 
is related to the choice of the asymptotic velocity. In the present case, = 0.05 a.u.. 
This value is in fact too high and leads to a strong confinement and thereby to a poor 
description of the shrinking of the bound states during the interaction with the pulse. 
When i?oo = 0.03 a.u., the discrepancies disappear and both curves are in good agreement 
at low electron energies. However, if the size of the basis stays equal to 1500, the result 
obtained by solving the scaled TDSE exhibits some unphysical oscillations at high energy. 
In that case, the confinement is weaker requiring an increase of the size of the basis to 
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Figure 9: (Color online) Electron energy spectrum resulting from the interaction of atomic hydrogen 
with a strong low frequency laser field. The cosine square pulse of peak intensity 10 Watt/cm 2 
and frequency oj = 0.114 a.u. has a total duration of 20 optical cycles. The blue dashed line is 
the result obtained by solving the unsealed TDSE with a basis of 800 Coulomb Sturmians per 
electron angular momentum and a dilation parameter n = 0.4 a.u.. The full red line is the result 
obtained by solving the scaled TDSE with a basis of 600 Coulomb Sturmian functions per electron 
angular momentum and the same dilation parameter. The scaled wave packet is propagated until 
i cll d = 20000 a.u. of time where it is stationary. 



cover a more extended region of space. An alternative could be to decrease the value of 
the dilation parameter. However, this dilation parameter fixes the spatial resolution in the 
whole space covered by the basis. It is therefore clear that a multi-resolution technique in 
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which the resolution is increased locally and gradually around the origin is more appropriate. 

In Fig. 9, we consider the case of atomic hydrogen interacting with a 20 optical cycle 
cosine square pulse of peak intensity 10 14 Watt/cm 2 and frequency 0.114 a.u.. The blue 
dashed line is the electron spectrum obtained by solving the unsealed TDSE with a basis of 
800 Coulomb Sturmian functions with a dilation parameter k equal to 0.4 a.u.. This result 
is in perfect agreement with the one obtained by Grum-Grzhimailo et al. [47] (see Fig. 4 
of that reference). The full red line is the electron spectrum obtained by solving the scaled 
TDSE with a basis of 600 Coulomb Sturmian functions with the same dilation parameter. 
Note that we have to propagate the scaled wave packet until t en d = 20000 a.u. i.e. during a 
long time after the end of the pulse. In fact, we have to wait until the wave packet becomes 
stationary before calculating the energy spectrum. In this context, it is therefore crucial to 
extract the scaled bound states. We clearly see that the results obtained by time scaling 
the coordinates are in perfect agreement with those obtained without scaling except at very 
small electron energies (first peak) where we observe a tiny difference which as before, is 
due to a slightly inaccurate description of the shrinking of the scaled bound states. In the 
present case, the scaling is switched on right at the beginning of the interaction and we use 
a small value of the asymptotic velocity, Roo = 0.001 a.u., because of the long duration of 
the pulse. 



V. CONCLUSIONS AND PERSPECTIVES 

In this contribution, we develop an ab initio approach to solve numerically the 
time- dependent Schrodinger equation that governs the ionization dynamics of atoms 
and molecules interacting with pulsed radiation fields. The approach is based on the 
combination of the time scaled coordinate method with an efficient time propagator. The 
key points of the time scaled coordinate method is a time-dependent scaling of the electron 
radial coordinate together with a phase transformation of the total wave packet of the 
system. This method presents the following advantages: (i) the fast oscillations resulting 
from the rapidly growing phase gradients are removed from the total wave packet thanks 
to the phase transformation, (ii) the scaled wave packet stays spatially confined while 
reaching a stationary state a sufficiently long time after the interaction with the pulse 
and (iii) the electron energy distribution is proportional to the modulus square of the 
scaled wave packet once it becomes stationary. This method has however an important 
drawback: it introduces different length scales in the problem. In particular, it leads to a 
shrinking of the scaled bound states. In principle, such an effect can be described properly 
by using a denser grid or a much bigger basis of C 2 functions. However, this inevitably 
increases the stiffness of the system of first order differential equations to solve for the 
time propagation of the scaled wave packet. Here, we show one efficient way of treating 
these problems. First, it is important to subtract the scaled bound states from the total 
wave packet after the end of the pulse, once the harmonic potential that confines the 
wave packet has disappeared. Second, we introduce a new high order time propagator of 
predictor-corrector type that so far revealed to be very efficient in handling the stiffness 
problem. The predictor is the fifth order explicit method of Fatunla and the corrector, a 
fully implicit Radau method of order seven. Despite the implicit character of the corrector, 
we show that all the calculations reduce to simple matrix-vector products that allow a high 
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level parallelization of the computer codes. Finally, our calculations suggest that an elegant 
way to further improve the efficiency of our method is the use of multi-resolution techniques. 

At this stage, the method has been tested in the case of the interaction of a pulsed 
radiation field with a one-dimensional model atom described by a Gaussian potential and 
with atomic hydrogen. Electron energy spectra have been calculated in rather demanding 
physical situations. In all cases, the new approach give very accurate results, particularly 
for high photolectron energies, at the expense of less computer resources when compared to 
the usual grid or spectral methods without scaling. 
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